started: AL13Jul2019
last updated: AL14Jul2019
Input sequencing data: 1,850 vars x 715 samples ( 519BC = 258UBC + 260CBC and 197NFE)
Input eigenvectors: 3,219 samples = 715 wecare-nfe + 2,504 kgen
Sys.time()
## [1] "2019-07-14 21:21:34 BST"
rm(list=ls())
graphics.off()
library(knitr)
## Warning: package 'knitr' was built under R version 3.5.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
base_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s14_joined_ampliseq_1kg_PCs/s02_not_rare_variants_not_in_LD_293"
opts_knit$set(root.dir = base_folder)
options(stringsAsFactors = F)
options(warnPartialMatchArgs = T,
warnPartialMatchAttr = T,
warnPartialMatchDollar = T)
#options(error = browser()) # Type Q or c to exit, drop browser level
# https://support.rstudio.com/hc/en-us/articles/200713843?version=1.1.456&mode=desktop
# https://stackoverflow.com/questions/13052522/how-to-leave-the-r-browser-mode-in-the-console-window/13052588
source_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s12_check_BRCA1_BRCA2_PALB2"
load(paste(source_folder, "s03_exclude_BRCA1_BCRA2_PALB2_carriers.RData", sep="/"))
base_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s14_joined_ampliseq_1kg_PCs/s02_not_rare_variants_not_in_LD_293"
eigenvectors_file <- paste(base_folder, "ampliseq_1kg_293_3219_100PCs.eigenvec", sep="/")
eigenvectors.df <- read.table(eigenvectors_file, header=T, sep="\t",quote="")
eigenvalues_file <- paste(base_folder, "ampliseq_1kg_293_3219_100PCs.eigenval", sep="/")
eigenvalues.df <- read.table(eigenvalues_file, header=F, sep="\t",quote="")
kg_phenotypes_file <- paste(base_folder, "integrated_call_samples_v3.20130502.ALL.panel", sep="/")
kg_phenotypes.df <- read.table(kg_phenotypes_file, header=T)
rm(source_folder, eigenvectors_file, eigenvalues_file, kg_phenotypes_file, genotypes.mx, variants.df)
ls()
## [1] "base_folder" "eigenvalues.df" "eigenvectors.df" "kg_phenotypes.df" "phenotypes.df"
dim(eigenvectors.df)
## [1] 3219 102
dim(eigenvalues.df)
## [1] 100 1
dim(kg_phenotypes.df)
## [1] 2504 4
dim(phenotypes.df)
## [1] 715 24
rownames(eigenvectors.df) <- eigenvectors.df$FID
eigenvectors.df <- eigenvectors.df[,-1]
2504+715
## [1] 3219
plot(eigenvalues.df$V1, type="b", ylab="Variance",
main="Top 100 eigenvectors")
plot(eigenvalues.df$V1[1:10], type="b", ylab="Variance",
main="Top 10 eigenvectors")
rm(eigenvalues.df)
# Make a table with IDs of overlapping NFE
eigenvectors.df[c(3022,3023),c("IID","PC1")]
## IID PC1
## 9_S346_L008 9_S346_L008 -0.0130949
## 2:HG00097 2:HG00097 -0.0085805
nfe_pca <- eigenvectors.df$IID[3023:3219]
nfe_ampliseq <- sub("2:","",nfe_pca)
nfe.df <- data.frame(nfe_ampliseq, nfe_pca)
# Remove overlapping NFE from eigenvectors data
selected_samples <- ! eigenvectors.df$IID %in% nfe.df$nfe_pca
sum(selected_samples)
## [1] 3022
518+2504
## [1] 3022
eigenvectors_ampliseq_kgen.df <- eigenvectors.df[selected_samples,1:6]
"sample" -> colnames(eigenvectors_ampliseq_kgen.df)[1]
# Prepare phenotypes for ampliseq-kgen data
phenotypes.df[c(518:519),c(1,2)]
## long_ids illumina_id
## 9_S346_L008 9_S346_L008 S346
## HG00097 HG00097 <NA>
phenotypes_ampliseq.df <- phenotypes.df[1:518,c("long_ids","cc")]
table(phenotypes_ampliseq.df$cc)
##
## 0 1
## 258 260
"WECARE" -> phenotypes_ampliseq.df$cc[phenotypes_ampliseq.df$cc==1]
"WECARE" -> phenotypes_ampliseq.df$cc[phenotypes_ampliseq.df$cc==0]
table(phenotypes_ampliseq.df$cc)
##
## WECARE
## 518
c("sample","group") -> colnames(phenotypes_ampliseq.df)
phenotypes_kgen.df <- kg_phenotypes.df[,c("sample","super_pop")]
c("sample","group") -> colnames(phenotypes_kgen.df)
phenotypes_ampliseq_kgen.df <- rbind(phenotypes_kgen.df,phenotypes_ampliseq.df)
table(phenotypes_ampliseq_kgen.df$group)
##
## AFR AMR EAS EUR SAS WECARE
## 661 347 504 503 489 518
# Add eigenvectors to phenotypes
eigenphen_ampliseq_kgen.df <- full_join(
phenotypes_ampliseq_kgen.df, eigenvectors_ampliseq_kgen.df, by="sample")
dim(eigenphen_ampliseq_kgen.df)
## [1] 3022 7
head(eigenphen_ampliseq_kgen.df)
## sample group PC1 PC2 PC3 PC4 PC5
## 1 HG00096 EUR -0.01345630 -0.01863560 0.01014450 -0.018405500 0.04524020
## 2 HG00097 EUR -0.00802251 -0.00903830 0.01379440 -0.021585900 0.01081080
## 3 HG00099 EUR -0.00887617 -0.00759887 0.02047870 0.010014300 0.01332560
## 4 HG00100 EUR -0.00823478 -0.01713840 0.01073160 0.006955610 0.02205440
## 5 HG00101 EUR -0.01322770 -0.01052920 0.01053310 -0.000193446 0.02723150
## 6 HG00102 EUR -0.00801074 -0.00603827 -0.00808896 -0.010795500 -0.00336974
tail(eigenphen_ampliseq_kgen.df)
## sample group PC1 PC2 PC3 PC4 PC5
## 3017 95_S517_L008 WECARE -0.0131126 0.000874616 0.00183998 -0.01053910 0.03642830
## 3018 96_S236_L007 WECARE -0.0104256 -0.009713690 0.01412110 -0.00247655 0.01483850
## 3019 97_S509_L008 WECARE -0.0117041 -0.016928600 0.02458060 -0.00367240 0.00259899
## 3020 98_S335_L008 WECARE -0.0117774 -0.011305600 0.00769248 -0.01457810 -0.00844976
## 3021 99_S418_L008 WECARE -0.0136038 -0.017416600 0.00380755 -0.00196730 0.00329286
## 3022 9_S346_L008 WECARE -0.0130949 -0.009817580 -0.00462172 -0.02054780 -0.01465900
# Clean-up
rm(nfe_pca, nfe_ampliseq, selected_samples, eigenvectors_ampliseq_kgen.df, phenotypes_ampliseq.df,
phenotypes_ampliseq_kgen.df)
# Prepare vector fr colour scale
myColours <- c("EUR"="BLUE", "AFR"="YELLOW", "AMR"="GREEN",
"SAS"="GREY", "EAS"="PINK",
"WECARE"="RED")
myColourScale <- scale_colour_manual(values=myColours)
# Static plot
ggplot(eigenphen_ampliseq_kgen.df, aes(PC1,PC2)) +
geom_point(aes(col=group)) +
labs(title="293 not rare variants not in LD", x="PC1", y="PC2") +
myColourScale
# Interactive plot
plotly_group <- factor(eigenphen_ampliseq_kgen.df$group,
levels=c("AFR","AMR","EAS","SAS","EUR","WECARE"))
g <- ggplot(eigenphen_ampliseq_kgen.df, aes(PC1,PC2)) +
geom_point(aes(col=plotly_group, text=sample)) +
labs(title="293 not rare variants not in LD", x="PC1", y="PC2") +
theme(legend.title=element_blank()) + # To suppress the legend title
myColourScale # otherwise it would be "plotly_group"
## Warning: Ignoring unknown aesthetics: text
ggplotly(g, tooltip="text") # By default the tooltip would also show coordinates
## Warning in dev_fun(file = tempfile(), width = width %||% 640, height = height %||% : partial argument match of 'file' to 'filename'
# Clean-up
rm(myColours, myColourScale, g)
ls()
## [1] "base_folder" "eigenphen_ampliseq_kgen.df" "eigenvectors.df" "kg_phenotypes.df" "nfe.df" "phenotypes_kgen.df" "phenotypes.df" "plotly_group"
sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS 10.14.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] plotly_4.9.0 ggplot2_3.2.0 dplyr_0.8.1 knitr_1.23
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 later_0.8.0 pillar_1.4.1 compiler_3.5.1 tools_3.5.1 digest_0.6.19 jsonlite_1.6 evaluate_0.14 tibble_2.1.3 gtable_0.3.0 viridisLite_0.3.0 pkgconfig_2.0.2 rlang_0.3.4 shiny_1.3.2 crosstalk_1.0.0 yaml_2.2.0 xfun_0.7 withr_2.1.2 stringr_1.4.0 httr_1.4.0 htmlwidgets_1.3 grid_3.5.1 tidyselect_0.2.5 glue_1.3.1 data.table_1.12.2 R6_2.4.0 rmarkdown_1.13 purrr_0.3.2 tidyr_0.8.3 magrittr_1.5 promises_1.0.1 scales_1.0.0 htmltools_0.3.6 assertthat_0.2.1 xtable_1.8-4 mime_0.7 colorspace_1.4-1 httpuv_1.5.1 labeling_0.3 stringi_1.4.3 lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4
Sys.time()
## [1] "2019-07-14 21:21:36 BST"